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CN \ Abstract 

We consider a class of finite Markov moment problems with arbitrary 
, number of positive and negative branches. We show criteria for the exis- 

fence and uniqueness of solutions, and we characterize in detail the non- 
unique solution families. Moreover, we present a constructive algorithm 
, to solve the moment problems numerically and prove that the algorithm 

computes the right solution. 



1 Introduction 



^ . We aim at inverting a moment system often associated with the prestigious 

■ name of Markov. The original form of the problem is the following. Given a 
T-H , finite set of moments mk for k = \,. . . ,K, find a bounded measurable density 

function / satisfying 

nik^ [ x''-^f{x)dx, 0</<l, k = l,...,K. (1) 

oo . 

, Condition for the existence of solutions f{x) to this problem is classical [1, 2]. 

^ ' In general solutions are not unique, unless more conditions are given, e.g. based 

, on entropy minimization [3, 4] or L°°-minimization [19, 18]. A typical result 

■ is that the unique solution for even K is piecewise constant, taking values in 
^ , {0, 1}. More precisely, it K — 2n then / is of the form 

_ _ I 

n 

where Xii^) is the characteristic function for the interval / and 

yi < Xi < 1/2 < X2 < ■ ■ ■ < yn < Xn- (3) 

See Theorem 3 below in Section 4 and consult e.g. [5, 8, 17, 23, 25] for general 
background on moment problems. 
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A reduced form of the finite moment problem is to search for sohitions to (1) 
which are precisely of the form (2, 3). One then obtains an algebraic problem 
for the branch values, 



1 

i'nk = T^x'^ -Vj, k = l,...,K = 2n. 



(4) 



Finding {xj} and {jjj} from {m.k} is an ill-conditioned problem when the branch 
values of the solution come close to each other; the Jacobian of the problem 
is a Vandermonde matrix and iterative numerical resolution routines require 
extremely good starting guesses when the matrix degenerates. For less than four 
moments a direct method based on solving polynomial equations was presented 
in [21]. Routines based on the Simplex algorithm have been proposed in [19]. 
Another algorithm was presented by Koborov, Sklyar and Fardigola in [16, 24] in 
the slightly modified setting where / takes values in {—1, 1} instead of {0, 1}. It 
consists of solving a sequence of high degree polynomial equations, constructed 
through a rather intricate process with unclear stability properties. In [14] we 
showed that this algorithm can be drastically simplified and adapted to (4). 
Later, in [15], we also gave a direct proof that the simplified algorithm indeed 
computes the correct solution, relying on the classical Newton's identities and 
Toeplitz matrix theory. 

The moment problem has many applications in for instance probability and 
statistics [10, 7], but also in areas like wave modulation [6, 22] and "shape 
from moments" inverse problems [11]. Our own motivation comes from a quite 
different field, namely multiphase geometrical optics [3, 4, 12, 13, 14, 21]. In 
this application one needs to solve a system of nonlinear hyperbolic conservation 
laws. To evaluate the flux function in the partial differential equations (PDEs) 
a system like (4) must be solved. In a finite difference method this means that 
the system miist be inverted once for every point in the computational grid, 
repeatedly in every timestep. It is thus important that the inversion can be done 
fast and accurately; this difficulty has been a bottleneck in computations. In [14] 
we used the simplified algorithm mentioned above for numerical implementation 
inside a shock-capturing finite difference solver. It is our aim here to develop 
better algorithms and understanding to open the way for the processing of 
intricate wave-fields with large K, and thus complement the seminal paper [4] 
where the multiphase geometrical optics PDEs were first proposed. 

In this paper we are concerned with a generalization of (4). In the geometri- 
cal optics application, the number of moments K is typically not even and one 
can have a variable number of positive (xk) and negative (yk) branches. We 
thus consider the following problem 



where Ux + riy = K but where rix and Uy arc not necessarily equal. We study 
existence and uniqueness of solutions to this problem (Theorem 2). In partic- 
ular we are interested in how and when uniqueness is lost. For these cases we 
characterize the family of solutions that exists. The reason is to understand 
what happens numerically close to degenerate solutions, which is an important 





(5) 
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feature in the application we have in mind: In the exact sohition to the muhi- 
phase geometrical optics PDEs the moment problem is typically degenerate for 
large domains; the numerical approximation is almost degenerate. 

We also give constructive algorithms to solve (5) and prove that they gener- 
ate the right solution (Theorem 1). In a future paper we will study the numer- 
ical stability of these algorithms. Experimentally we note, for instance, that to 
compute the next moment. Algorithm 3 is much more stable than Algorithm 1. 
The difHculty lies in understanding perturbations around degenerate solutions, 
which is where the algorithms are most unstable. For this the insights of this 
paper will be of importance. 

Remark 1 The problem (5) can be cast in the form of (1) if one demands that 
the density function f{x) is of the form 

f{x) = J2 sgnfe) mx) - H{x - \xM - E ^^"^^y^) [^(^) - - l^'^ l)] ' 

(6) 

and we rescale the moments ruk km^. For the case rix = riy = n and K = 2n 
with interlaced branch values (3) this reduces to (2). 

This paper is organized as follows. In Section 2 we present the algorithms for 
solving (5). Notation and various ways of describing a solution is subsequently 
introduced in Section 3. Next we derive conditions for existence and uniqueness 
of solutions in Section 4 and also discuss various properties of the solution, in 
particular when it is not unique. A theorem proving the correctness of the 
algorithms is proved in Section 5. Finally, in Section 6, we give additional 
properties of the elements of our algorithms, and use these to relate our results 
back to the classical Markov theory. 



2 Algorithms 

In this section we detail the algorithms that we propose for solving (5). The 
solution that we obtain is what we call the minimal degree solution, meaning 
that when the solution is not unique as many branch values as possible are zero. 
See Section 4 for a precise definition. The algorithms goes as follows; they may 
fail in case there is no solution to (5). 

Algorithm 1 (Computing {xj} and {%}) 

1. Construct the sequence {ofe} as follows. Set oq = 1 and Ofe = for k < 0. 
For 1 < k < K , let the elements be given as the solution to 
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2. Construct the matrix Ai e ]R"i«^"» as 



( O-ny ^ny-\ ■■■ (^ny-na;+l\ 

luy + l 0,ny ■■■ Ctn„-n^+2 



y^Uy+n^— 1 «nj,+n^-2 ■•• dny J 

Compute the rank of Ai . Let fix = rank A\ and hy = Uy — Ux + fix- 
3. Construct the matrices Ao,Ai e ]R"»><"« as 



Ao 



Ai = 



( any+\ afiy 

O'ny+2 0-ny+\ 

I any ttny-l 

"■n„ + l n.ny 

\0'hy+h:c-l 0,ny+h^-2 



0'ny-h^+2\ 
afiy-nx+S 

afiy + l J 

dfiy-n^+l \ 

dhy-n^ + l + l 

any ) 



4- Solve the generalized eigenvalue problem 

Aov = xAiv, (8) 

to get the {xj} values of the minimal degree solution to (5). 

5. To compute the {yj} values, the same process is used with ruk replaced by 
—ruk and the roles of Ux and Uy interchanged. 

An alternative to Algorithm 1 is as follows: 

Algorithm 2 (Computing {xj} and {yj}) 

1. Construct the matrices Aq and A\ as in steps 1-3 in Algorithm 1. 

2. Denote the first column vector in Aq by do by and solve 

Aic' = -ao, c' = {Ci,C2,...,Cnj'^. (9) 

3. Construct the polynomial 



P{z) = Cn^ + Cii^_iZ + 



The roots of P{z) are the {xj} values of the minimal degree solution to 
(5) (possibly together with some zeros). 

4. To compute the {yj} values, the same process is used with ruk replaced by 
—rUk and the roles of Ux and Uy interchanged. 



Remark 2 We note that the values of Ok in the definition ( 7) are independent 
of K, since the system matrix is triangular. We therefore consider the sequence 
without reference to K in any other respect than the fact that we are only able 
to com.pute elements with k < K tuhen we are given K moments. The largest 
index of the ak-sequence a,ppea,ring in the matrix Ai is Uy + Ux — i < K . In the 
matrices Aq,Ai it is riy + = Uy — + 2nx < riy +nx = K. Hence all three 
matrices can he constructed from the first K moments. Some properties of the 
A\ matrix are detailed in Section 6. 

Sometimes one is not interested in finding the individual {xj} and {j/j} 
branch values but just wants the higher moments, defined as 

mk = Y.^'i-Y.yl (10) 

but now for k > K, given a solution {xj}U{yj} to (5). (That this is well-defined 
is shown later in Theorem 2.) For this case there is another algorithm, which 
has empirically proven to be more stable than first computing {xj} and {yj} 
from Algorithm 1 or 2, and then entering the values into (10). We stress that 
this is precisely what is needed in order to compute /C-multivalued solutions of 
the inviscid Burger's equation in geometrical optics, following the ideas of [4]. 

Algorithm 3 (Computing niK+i) 

1. Construct the A\ matrix as in steps 1-2 of Algorithm 1. 

2. Let 

ClO = i<lny+l, ClrLy+2, ■ ■ ■ , dn.y+n^)'^ G M"=^ , 

and let c = {ci,C2, . . . , be one solution to 

Aic=-ao. (11) 

3. The next moment is given by 

K 

ruK+i = -{K + l)^CjaK+i-o - ^m3<^K+i-o- 

3=1 3=1 

We recall that Algorithm 1 has been shown to be numerically efficient in the 
paper [14]. The justification of these algorithms is given in Section 5 where we 

show the following theorem: 

Theorem 1 If a solution to (5) exists then: 

(i) In Algorithm 1, the matrix Ax is non-singular. The generalized eigenvalue 
problem in (8) is well-defined and the generalized eigenvalues (counting 
algebraic multiplicity) are the {xj} -values of the minimal degree solution 
to (5) plus fix — Drain zeros. (See (19) for the definition of D^in-) 

(ii) In Algorithm 2, c' is well defined, 

P{z) = det{zl - A^^Ao) (12) 

and the roots of P{z) are the {xj}-values of the minimal degree solution 
to (5) plus fix 
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(iii) In Algorithm 3, the computed moment satisfies 



i=i j=i 

for all solutions {xj} U {yj} to (5). 

We postpone the proof of Theorem 1 to Section 5. We just note here that the 
last point in Algorithms 1 and 2 can easily be explained by the symmetry of the 
problem. Indeed, the negative of (5) 

-mfe = ^y)-24' k = l,...,K, 

is of the same form as (5) itself, with the roles of Hx, {xj} and Uy, {yj} inter- 
changed. 



3 Preliminaries 

We will use three different ways of describing the solution to (5). First we 
have a set of numbers {xj}"^^ and solving (5). We call those numbers 

branch values. Second, we have a pair of polynomials {p, q) of degrees at most 
rix and Uy respectively in the z variable. Third, wc have a pair of coefficient 
vectors c = (co,...,c„J^ e and d = (do, • • • , c?nj^ e These 

three representations are related as 

p{z) = (1 - xiz) •••(!- Xn^z) = Co + CiZ + • • • + c„^_i2;"=^-i + Cn.z"", (13) 

and 

q{z) = (1 - yiz) ■••(!- y^^z) = do + diz + ■ ■ ■ + dn^-iz^'y'^ + dn^z'^K (14) 

It is clear that there is a one-to-one correspondence between these ways of 
describing the solution, if wc disregard the ambiguity in the ordering of the 
numbers {xj} and {yj}- Generally, we will use the notation Deg(p) to denote the 
degree of a polynomial p, and, for a given coefficient vector c, we systematically 
write Pc to denote the corresponding polynomial (13). 

Definition 1 We call the pair of polynomials {p,q) a (polynomial) solution to 
(5) if 

1. The degrees of p and q are at most and Uy, 

Deg(p) < Ux, Deg(g) < Uy, (15) 

2. They are normalized to one at the origin, 

p(0) = g(0) = 1, (16) 
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3. Their roots {xj} and {yj} satisfy 

Deg(p) Deg(q) 

rnk= ^t- E k = l,...,K. (17) 

j=l 3 = 1 

We note that the roots cannot be zero because of (16). 

Next: 

Definition 2 A pair of vectors 

c = (Co, . . . , c„J^ e and d = {do, • • • , G 

is said to be a (coefficient) solution to (5) if the corresponding pair {Pc,Pd) 
(13)-(14) realizes a polynomial solution to (5). 

The number of branch values are always Ux and Uy respectively. Some of 
them may be zero, and they do not need to be distinct. The number of non-zero 
branch values are Deg(p) and Deg(q') respectively. The degree of a solution can 
then also be defined. 

Definition 3 The degree of a solution to (5) is the number of non-zero xj- 
values. This number is equivalent to Deg(p). 

Given any polynomial pair satisfying (16), we say that it generates the mo- 
ment sequence {nik} if rrife is given by (17) for all k. In turn, each sequence 
of moments {rrik} generates the corresponding {uk} sequence through (7). We 
define the big matrix 



We let the columns of A be denoted ao, . . . , and we note that 

'I I W I Wl \ 

a-o ■■■ I = I ^0 an\ = lao Ai I, (18) 

Hence, Aq and Ai constitutes the first and last columns of A respectively. 
When ao € range Ai and ao 0, let 

-Dmin = argmin^>o oq G span{ai, . . . , aj}, (19) 

and set Dmin = if ao = 0. Moreover, define 

Ux — rank A\. (20) 



4 Existence and uniqueness of solutions 



In this section we prove results on the existence and uniqueness of solutions to 
(5). We aim at establishing the following theorem: 

Theorem 2 

(i) There exists a solution to (5) if and only if 

ao e range(Ai). (21) 

(ii) If d is the degree of a solution to (5), then .Dmin ^d < -Dmax- 

(iii) When (21) holds, there is a unique solution {p*,q*) of minimal degree 
-Dmin- For this solution, xj ^ j/j for all indices i,j representing non-zero 
branch values. Moreover, Deg{q*) < riy — + rank Ai with equality if 
-Dmin < rank Ai . 

(iv) When (21) holds, a polynomial pair {p,q) is a solution if and only if 
p = p*r and q = q*r where r{z) is a polynomial satisfying r(0) = 1 
and Deg(r) < £>max - -Dmin- 

(v) The minimal degree solution is the only solution to (5) if and only if the 
matrix Ai is non-singular. 

(vi) Let {xj} and {yj} be a solution to (5). Then the higher moments defined 
in (10) are well-defined. 

Let us proceed with several remarks: 

Remark 3 In particular it follows from (i) that there exists a solution as soon 
as the matrix Ai is non-singular. 

Remark 4 Since (5) is a system of polynomial equations of degree K , one could 
expect there to be a finite number of solutions, typically K solutions. How- 
ever, because of the special structure of the equations there is either one unique 
solution (when Ai is non-singular) or inifintely many solutions (when Ai is 
singular). 

Remark 5 The form {p*r,q*r) of solutions can also be stated as follows: All 
solutions have a core set of values {xj}, j = 1. . . . ,Deg(p*) = -Dmin ond {yi}, 
i = 1,. . . ,Deg{q*) corresponding to non-zero branch values of the minimal degree 
solution, where xj ^ yi for all those One can then add an optional set of 
non-zero branch values {a;£)„i„+j}, and {2/Deg(g*)+j}; for j = 1,. ..,D 

max -^min 

such that XD^i^+j = yBes{q*)-\-j- 

To prove this theorem we first establish some utility results in the next sub- 
section. We then derive different ways of characterizing the solution in Section 
4.2, which are subsequently used to prove Theorem 2 in Section 4.3. 
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4.1 Utility results 

We start with a useful lemma on Taylor coefRcients for a product of functions: 

Lemma 1 Suppose f, g and h are analytic functions in a neighborhood of zero 
satsifying f{z) = g{z)h{z). Let f have the Taylor expansion 

oo 

f{z) = Y.fkz', 

fc=0 

and let {gk}, and {hk} be the corresponding coefficients for g{x) and h{x) re- 
spectively. Then 

k 

fk^J293>^k-j. (22) 
Proof: Since the functions are analytic the coefRcients are given as 



j=0 



where Cjk = kl/j\{k — j)! arc the binomial coefficients. But g^^\0) = j\gj and 
/i('=-j)(0) = {k-j)\hk-j and therefore (22) follows. □ 

Remark 6 The sum (22) is in fact precisely an elementwise description of 
multiplication of a lower triangular k x k Toeplitz matrix by a vector. In the 
notation of [15], it would read f = T{g)h. 

As was already known by Markov, the exponential transform of the moment 

sequence plays an important role in the analysis of these problems, sec e.g. [1, 2]. 
We show here that {a^} is a version of the exponential transform of {to^}. 

Lemma 2 Suppose {m^} is generated by the polynomials p{z) and q{z) and 
{ak} is generated by {ruk}. Let m{z) be defined as 

m{z) = miz + ^m2Z^ + ^™3-2^ H • (23) 



Then if (16) holds, 

^ 9(f) =ao + ai^ + a2^' + --- , (24) 



p{z) 

written as its Taylor expansion around z = 0. 

Proof: Let us first show that m{z) is a well-defined analytic function at zero. 
We have 



^ ^ fc fc 

k=0 l^j=l ~~k l^k=0 ^j=l ~~k~ 
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The last step is allowed when \z\ < l/maxij(|xj|, which is true for small 
enough z since p(0) ^ 0. This also shows that the function is analytic at zero. 
Moreover, 

Finally, setting a{z) := exp(m(2:)) and differentiating gives 

za'{z) = zm'{z)a{z), 

where all three functions arc analytic at zero. Let a{z) have the Taylor co- 
efficients {cik}- Then za'{z) — aiz + ^a^z^ + Za-^z^ ■ ■ ■ and clearly zm!{z) = 
miz + m2Z^ + • • • . By Lemma 1, for k > 1, 

k 

j=i 

Since fio = q{0)/p{0) = 1, wc sec that Ofc and dk satisfy the same non-singular 
linear system of equations (7), and therefore Uk = dk, showing (24). □ 
We now have the following basic characterization of a solution. 

Lemma 3 Suppose p{z) andq{z) are two polynomials satisfying (15, 16). They 
form, a polynomial solution to ( 5) if and only if their quotient has the Taylor 
expansion around z = 

44=«0 + ai^ + --- + ax^^ + O(^^+i), (25) 

p{z) 

where {a^} is generated by {wfe}. Moreover, if {p,q) is a solution then {p,q) 
is also a solution if and only if the pair satisfies (15, 16) and p/q = p/q where 
these fractions are defined. 

Proof: Let {fhk] be generated by p and q and suppose (25) holds. Then, as in 
the of proof of Lemma 2 for 1 <k < K 

k 

kok = y^^rhjak-j- 
Since {ruk} satisfy the linear system (7), we have after subtraction, 

rUn - fhn = -'^(nik - rhk)an-k, mi = rhi, 

fe=i 

for n = 2,. . . ,K. By induction ifik = ruk for 1 < k < K, showing that {p, q) 
solves (5). On the other hand, if {p, q) is a solution, then (25) must hold by (24) 
in Lemma 2. 

For the last statement, the "if" part is obvious since both pairs then satisfy 
(25). To show the "only if" part, suppose both {p,q) and {p,q) are solutions. 
By definition they satisfy (15, 16), and by (25), 

q{z) q{z) _ q{z)p{z) - p{z)q{z) _ ^,k+i^ 



p{z) p{z) p{z)p{z) 



= 0{z' 
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Since p(0)p(0) = 1 we must have that {q{z)p{z) — p{z)q{z))/ z^^^ is bounded 
as z — > 0. But since the degree of qp — pq is at most K = nx + riy this is only 
possible if it is identically zero. Hence q{z)p{z) = p{z)q{z) which concludes the 
proof. □ 

4.2 Characterization of the solution 

In this section we show three Propositions that characterize solutions to (5) 
in terms of polynomials, coefficient vectors and the column vectors of the A- 
matrix in (18). Wc start by expressing the uniqueness properties of the solution 
in terms of its polyomial representation. 

Proposition 1 Suppose the pairs {p, q) and {p, q) are both polynomial solutions 

to (5). Then, 

(i) Deg(p) - Deg(g) = Dcg(p) - Deg(g). 

(ii) //Deg(p) < Deg{p), and if there is no polynomial r{z) such that p = pr, 
then there is another solution {p,q) with Deg(p) < Deg(p). In particular, 
if Deg(p) = Deg{p) but p^p, there is such a lower degree solution. 

(iii) If Deg{p) < Deg{p), any polynomial pair {pr,qr) is a solution if r{z) 

is a polynomial satisfying r(0) = 1 and Dcg(r) < Dcg(j)) — Dcg(p). In 
particular, ifDeg{p) < m < Deg(p) there is a solution {p, q) with Deg(p) = 
m. 

Proof: 

(i) The statement follows directly from Lemma 3, since qp = pq implies that 

Deg(g) + Deg(p) = Deg(p) + Deg(g'). 

(ii) We let 

P{z) = rp{z)p{z) + Sp{z), q{z) = rq{z)q{z) + Sq{z), 

be the unique polynomial decomposition of {j),q) such that rp,rq,Sp,Sq 
are polynomials, Deg(sp) < Deg(p) and Deg(sg) < Deg(g). Since pq = pq 
by Lemma 3, we get 

Pqirq -rp) = qsp - psg. 

Unless Tg = rp the degree of the left hand side is at least Deg(p) + Deg(g), 
while the degree of the right hand side is at most 

max (Deg(g) + Deg(sp), Deg(p) + Deg(sg)) < Deg(g) + Deg{p). 

Hence, = rp and qSp = pSq. Since q,p ^ it follows that either Sp and 
Sq are both zero or both non-zero. Suppose Sp ^ and Sq ^ 0. Write 
Sp{z) = z"^''Sp{z) and Sq{z) = z"^'>Sq{z) where Sp(0) 7^ and Sg(0) 7^ 0. 
Since 

z"'''Sp{z)q{z) = z"'"'Sq{z)p{z) 
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and also q{0) = p{Q) = 1, the lowest degree term in the left and right hand 
side polynomials are z"^" and z™" respectively, and therefore nip = m,. 
Consequently, 

Sp{z)q{z) = Sq{z)p{z), 

and Sp(0) = 5^(0). We can then take p{z) = Sp{z)/sp{0) and q{z) = 
Sq{z)/sq{0). They satisfy 

p{z)q{z) = q{z)p{z), p{0) = q{0) = 1, 

while Deg(p) = Deg(sp) < Deg(sp) < Deg(p) and similarly Deg((J) < 
Deg{q) < Uy. Hence {p,q) is a polynomial solution by Lemma 3. It 
has degree strictly less than {p,q), which shows the first statement in 

(ii). If Deg(p) = Deg(p) and p ^ p then there is no r{z) satisfying the 
requirements, showing the second statement in (ii). 

(iii) We finally let r{z) be any polynomial with Deg(r) < Deg(p) — Deg(p) and 
r(0) = 1. We then set p = pr and q = qr. These polynomials trivially 
satisfy (16) and (25). Since Deg(p) = Deg(r) + Deg(p) < Deg{p) < 
and 

Deg(g) = Deg(r) + Deg(g) < Deg(p) - Deg(p) + Deg(g) = Deg(g) < %, 

they also satisfy (15) and thus are a polynomial solution by Lemma 3. In 
particular we can take r{z) of degree m. 

□ 

A solution to (5) can also be characterized in terms of the coefficient vectors. 
We have the following Proposition. 

Proposition 2 The pair c = (cq, . . . , c„^)-^ e M"»'+^ and d = {do, . . . , dn^y € 
M"s'+^ is a coefficient solution to (5) if and only if 

(i) Co = 1, 

(ii) c is in the null-space of A, 
(iii) 

min(fc,n^ ) 

dk= ^ CjUk-j, k = 0,...,ny. (26) 

3=0 

Proof: Suppose first that c is in the null-space of A, cq = 1 and {dk} is given by 
(26). Extend the coefficient sequences by setting = for fc > rix and dk = 
for k > Uy. Since c is in the null-space of A, we get '^j^QCjUk-j = when 
Hy + 1 < k < Ux + Ux = K, and in conclusion 

k 

dk^^Cjak-j, k = 0,...,K. (27) 

Upon noting that {cfcj^Q and {c?fe}^o ^^c Taylor coefficients of Pc and Pa, 
and since Pc{0) = Cq = 1, Pd{0) = c?o = aoco = 1, Lemma 1 shows that 

Prf(z) = Pc{z) [ao + aiz + --- + gkz'' + O {z''+')] , (28) 
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and by Lemma 3 wc have that (Pc, Pd) is a solution to (5). Conversely, if {Pc, Pd) 
is a solution, then cq = -Pc(O) = 1 and by Lemma 1 we get that (27) holds. For 
k = Hy . . . ,K this also implies that c is in the null-space of A. □ 

The final Proposition of this section relates the degree of the solution to the 
column vectors of A and the linear spaces they span. 

Proposition 3 Let Vj = spanjai, . . . , Oj} and Vj^ = spanjao, . . . , aj}. Set 

Vo = V^i=9. Then 

(i) There is a solution if and only if ao G Vn^ — Range(Ai). 

(ii) There is a solution of degree j >0 if and only if 

ao G Vj, and Uj e (29) 

(iii) When ao G V^^ then 

ao e Vd, = Vd, 

if and only if d> Umax- 

(iv) When ao G Vn^ the vectors 

ai, . . .,aD^i^, 

(when -Dmin > 0^ 

(when -Dmax < n^), are all linearly independent. Moreover, 

Oj G yD„i„, Vj = Vd^,„, j = -Dmin, • • • , -Dmax- 

Proof: 

(i) By Proposition 2 there exists a solution to (5) if and only if there is a 
coefficient vector c = (1, c')^ in the null-space of A, i.e. 

Ac = Aic + ao = 0. 

But such a vector c exists if and only if ao is in the range of Ai. This 
shows (i). 

(ii) Again by Proposition 2 there is a solution of degree j if and only if there 
is a vector c = (cq, ci, . . . , c^, 0, . . . , 0)-^ such that 

= Ac = coOq + ciOi + ■ • • + Cjaj, (30) 

with Cj ^ and cq = 1. For j = this is clearly equivalent to ao = 
ov ao & Vo = V^i- For j > the existence of -coefficients satisfying 
(30) is equivalent to the left condition in (29). Moreover, if aj ^ V^_^ = 
span{ao, . . . , Oj-i}, then we must have Cj = to satisfy (30), and c 
cannot represent a solution of degree j. On the other hand, if Cj = 
and aj = CqCo + ■ • ■ + Cj_^aj^i for some non-zero coefficients c'^,, then 
ao + c'{ai -!-••• + c'-_-^aj-i + aj = 0, with dj. = (1 + Co)cfe — cJ., represents 
a solution of degree j. This shows (ii). 
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(iii) The statement is obvious in case -Dmin = 0. If .Dmin > there are scalars 
such that 

ao = viai H h vo^i^ao^i^, (31) 

by (19). Hence, oq € Vd^^^ and since the Vj spaces are nested, Vj C T^+i, 

wc have ao G Vd for d > -Dmin- Moreover, the minimal property of Z^min 
ensures that VD^i„ ^ in (31), so that ao ^ Vd when d < -Dmin- 

(iv) To show that when -Dmin > the vectors Oi, . . . , ao^ir, are hnearly inde- 
pendent, we use (31) and note that Pc{z) with c = (1, — Ui, - . - , —VD^i^i 0) • • 
is a polynomial solution to (5). Suppose now that the there are non-zero 
coefficients c'j such that 

c'^ai H h c^„i„aD^,„ = 0. 

Then P^' with c' = (1, c'j — v\, . . . , Cd^i^ — f^D^in) 0; ■ ■ ■ 1 0)^ is another 
polynomial solution to (5). Moreover, by the minimality property of -Dmin 
we must have CD^i^—VD^i^ ^ and therefore Deg(Pc) = E)eg(Pc') = -Dmin- 
But by (ii) in Proposition 1 this implies that there is yet another solution 
Pc" of degree strictly less than -Dmin- Hence, there are coefficients c" such 
that 

ao + c'/ai H h c^a^ = 0, 

with d < -Dmin, contradicting (19). The vectors must therefore be linearly 
independent. 

Suppose D* > -Dmin is the highest degree of an existing solution. Since 
Pc(z) is a solution of degree -Dmin we get from (iii) in Proposition 1 that 
there are solutions of all intermediate degrees -Dmin, . - - , D* ■ Hence, from 
(ii), aj e Vj^-i for j = -Dmin, - - - , -D* and from (iii) G Vj-i for j = 
-Dmin+1, - - - , -D*. Noting that if a^+i e Vj then Vj ~ Vj+i we can conclude 
inductively that VD„i„ = - - - = Vd* and aj G Vd^j^ for j = -Dmin, . - . , -D*. 
We now have three different cases: 

1. If -D* = Ux then VD„,i„ = Vn^ and by (20) we get D* = rank Ai — 

•Dmin + D'max = dim V„^ --Dmin + -Dmax = dim VD„i„ - -Dmin + -Dmax = 

-Dmax since either ai, . . . , aD„i„ are linearly independent or -Dmin = 
and Vo^in = 0- This shows (iv) for D* =nx- 

2. If -D* < Ux and Dmin = then Vd^,^ = Fc,. = and 

Vn^ = span{a£).+i, . . . , a„J. (32) 
Suppose there are non-zero coefficients ak such that 
aD'+iao'+i H h Q!„^an» = 0, 

and let k* be the highest index of all non-zero coefficients, ak* ^ 0. 

Then a^* G V^,_-^ and there is a solution of degree k* by (ii), a 
contradiction to the definition of D*. Hence, the vectors in (32) 
must be linearly independent and 

D* =nx - dim Vn^ = Dmin + nx - rank Ai = -Dmax, 

showing (iv) for this case. 
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3. If D* < Tlx and -Dmin > we have 

Vn^ = span{ai, . . . ,aD^,„,aD.+i, . . . ,a„^}. (33) 

Suppose there are non-zero coefficients au such that 

aioi H h aD„i„aD„,„ H h Q;_D.+iaD.+i H h a„^a„^ = 0. 

Since ai, . . . ^a^^^^ are linearly independent at least one ak with 
k > D* must be non-zero. By the same argument as above in case 
two we then get a contradiction and the vectors in (33) must be 
linearly independent. Hence, 

D* = Dmin + nx- dim Vn^ = i'min + "a: " rank Ai = Dmax, 

showing this final case. 

□ 

4.3 Proof of Theorem 2 

To prove Theorem 2 we essentially have to combine the results from Propositions 
1 and 3. The statement (i) is given directly by (i) in the latter. For the remaining 

points wc have: 

(ii) From (ii) in Proposition 3 we see that ao € Vd and ad S V^_-^. It follows 
from (iii) in Proposition 3 that d > .Dmin. On the other hand, if Dmax < 
Ux and d > -Dmax it says that V,i_i = Vd-i- Hence, ad G Vd-i which 
contradicts the linear independence of a D„ax ) • • • ) *^n^ established in point 
(iv) of Proposition 3. 

(iii) We note that by (19) there are scalars V\,. . . , Vo^ir, such 

ao = t^itti H h VD^,„a£)„i„. (34) 

Hence, Oq G Vd,,,;,, and since ?;D„i„ 7^ 0, we also have aD„i„ G . . 
By (ii) in Proposition 3 there is thus a solution of degree I^min which we 
denote {p*,q*). Since ai, . . . , Od^j^ are linearly independent by (iii) in 
Proposition 3, the coefficients in (34) are unique and therefore also the 
-Dmin-degree solution is unique. Moreover, suppose that Xj = yi = x* ^ 
for some Then p* and q* would have a common factor (1 — zx*), and 
by Lemma 3 also p(z) := p*{z) /{I- zx*) and q{z) := q* {z)/{l- zx*) would 
be a solution. But this is impossible since Deg(p) < Deg(p*) = Dmin- By 
(iv), shown below, a solution is given by {p*r,q*r) where r(0) = 1 and 
Deg(r) = Dmax - Dmin- Heucc Uy > Deg(g*r) = Dcg(g*) + - rank Ai. 
Suppose finally that Dmin < rank Ai and that Deg(g*) < ny — rix + 
rank Ai. Let Deg(r) = Dmax + 1 — Dmin- Then {p*r,q*r) is still a 
solution by Lemma 3 since {p* ,q*) is a solution, Deg(p*r) = Dmax + 1 = 
Ux + Dmin + 1 — rank Ai < Ux and 

Deg(g*r) < Uy - Hx + rank Ai + Dmax + 1 - Dmin = Uy + 1. 

This contradicts (ii) and therefore Deg(g*) = %— n^j+rank Ai, concluding 
the proof of (iii). 
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(iv) Wc first note that there exists a solution of degree Dmax by Proposition 
3 since if Dmax > -Dmin we have Oq G ^B„a>.-i ^^'^ '^d^^^ € Vd„,„ = 
^iimax-i — ^i)max-i- Hcncc, (iii) in Proposition 1 shows that any poly- 
nomial pair of the stated type is a solution. On the other hand, if the 
polynomial solution is not of this type, then (ii) in Proposition 1 says 
there is a solution of degree strictly less than D^i^, contradicting (ii) 
above. 

(v) We suppose first that A\ is non-singular. Then rank Ai = rix so that 

Dmin = -Dinax ^nd thc uniqucncss is given by (iii) above. If, on the con- 
trary, Ai is singular then Z^max > -Dmin and since we can then pick in- 
finitely many polynomials r{z) in (iv), we have infinitely many solutions. 

(vi) This is a consequence of (iv). The solution can be represented by {p*r, q*r) 
for some polynomial r{z) with r(0) = 1. Let l/xj for j = 1, . . . , Dj^in and 
l/jjj for j = 1, . . . ,Dcg((j*) be the roots of p*(z) and q* (z) respectively. 
Let 1/zj for j = 1, . . . ,Deg(r) be the roots of r{z). Then 

Dn,i„ Dog(r) Deg(g*) Deg(r) D„i„ Deg(g*) 

"^fc = E 4 + E - E - E = E 4 - E yj' 
j=i j=i j=i j=i j=i j=i 

which is independent of r{z) and uniquely determined because {p*,q*) is 
unique. 

5 Proof of Theorem 1 

We can now use the results in Section 4 to prove Theorem 1. 

(i-ii) To show the statements about Algorithms 1 and 2 we consider the reduced 
problem 

fix fly 

"^fe = E*^'-E^i' k = l,...,K, (35) 

where hx = rank Ai < rix, fiy = ny — rix+fix l£ ny and K = fix+fiy < K . 
The moments ruk in the left hand side are the same as in (5). First, we 
consider the minimal solution {p*-,q*) of (5). By (iv) in Proposition 3 
we must have Deg(p*) = -Dmin < rank A\ = hx- Moreover, by (iii) in 
Theorem 2, 

Deg(g*) < % — Tlx + rank Ai =ny. 

It follows from Lemma 3 that {p*, q*) is also a solution to (35). Second, let 

(p* , q*) be thc minimal degree solution to (35). Then by (iv) in Theorem 
2 there is a polynomial r{z) with r(0) = 1 such that p* = p*r and q* = 
q*r. But then {p*,q*) is also a solution to (5) by Lemma 3. By the 
uniqueness of the minimal degree solution of (5) it follows that r = 1 and 
p* = p*,q* = q* ■ Suppose now that there is another polynomial r{z) with 
r(0) = 1, Deg(r) > such that {p*r,q*r) is a solution to (35). Then 
Deg(p*r) = Dmin + Dcg(r) <nx= rank Ai. Hence, Dmin < rank Ai and 
therefore by (iii) in Theorem 2 we have Deg(g*) = % — n^ + rank Ai = %. 
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Thus, Deg(g*r) > hy which is impossible if (p*r, q*r) is a solution. Hence, 
is the unique solution to (35) and therefore Ai is non-singular by 
(v) in Theorem 2. 

Since Ai is invertible, the generalized eigenvalue problem (8) and c' are 
well-defined. Moreover, we can construct A^^Aq. By (9), 



A^'Ao = 



( -ci 

-C2 







1 



which is a companion matrix. It is well-known that for those matrices 
the elements in the first column are the coefficients of its characteristic 
polynomial. This is shown as follows: let be the minor of Y := 
zl — AJ'^Aq, i.e. the determinant of the matrix obtained by removing row 
i and column j. Then, the determinant can be expanded by minors, for 
any j, 

det(y) = 5^(-l)'+%i,M,,-, V = {vij}. 

i=l 

Taking j = 1, we get Mj_i = det(diag(2;, ■ ■ ■ ,z, —1, . . . , —1)) with i — 1 
occurrences of —1, so that Mj^i = ^;"»'~*(— Therefore, 



det{V) 



(-l)2(ci + z)Mi,i + Er=2(-l)'+'ciM,,i 



= ClZ""'~ 



=2 '-^ 



which is exactly (12). This shows that the results of Algorithms 1 and 2 
are identical, since the generalized eigenvalues in (8) are exactly the roots 
oiP{z). 

It remains to show what the roots are. Let A = [ciq Ai] be the A-matrix 
related to (35). Clearly, c ~ (1, c'^)^ is in the null-space of A and hence 
Pc{z) is the unique solution to (35). But for z ^ 0, 



P{z 



-1Z+- 



■ Ciz 



z J 



= z^-{^ 

= 0"^Pe(l/^) 

= ^"==(1 - Xi/Z)[\ - X2lz) • • • (1 - XD^,JZ) 
"(^ - Xi){z -X2)---{Z- XD^tJ, 



= z 



which extends to 2; = by continuity. This concludes the proof of points 
(i) and (ii). 

(iii) Let {p, q) be a polynomial solution to (5) and c the corresponding coeffi- 
cient solution. Prom Lemma 2 we have 



g(z) = p(z)e™(^\ 
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where m(z) is defined in (23). For the (K + l)-th Taylor cocfRcient of the 
left and right hand side we have by Lemma 2 and Lemma 1, 



= '^CjaK+i-j => ttK+i = -^aK+i-jCj, (36) 

since the A'-th Taylor eoefBeient of q and p is zero for k > and k > Hy 
respectively. Finally, the last row of (7) extended to size K + 1 gives 

K 

rriK+i = [K + l)aK+i - ^mjaK+i-j- 
Together the last two equations show point (iii). 



6 Properties of Ai and Markov's Theorem 

We now look more in detail on the structure of the Ai matrix. In particular we 
look at the implications of AiR being positive definite. Then we get an explicit 
simplified formula for the matrix and our results also shed some light on the 
relationship of our results to the classical Markov theorem on the existence 
and uniqueness of solutions to the finite moment problem (1) discussed in the 
introduction. For this we need to define the matrix 



R = 



V 



and note that left (right) multiplication by R reverses the order of rows (columns) 
of a matrix. In our notation we can then formulate Markov's theorem as follows 

Theorem 3 (Markov) Suppose K = 2n is even and There is a 

unique piecewise continuous function f{x) satisfying 



rrik 



= k [ x''-'^f{x)dx, 0</<l, k = l,...,K, (37) 



if AiR is symmetric positive definite and the matrix 

( ao Ai 

T 

\aK+i clq 

is singular. This f is of the form in (2, 3). 



(38) 



Remark 7 The theorem, does not rule out other forms of f{x) a priori, and 
without the second condition in (38) such solutions are indeed possible. It only 
considers the case Ux = Uy, i.e. problem (4), and says nothing about the possi- 
bility of other solution types, e.g. when the {xj} and {yj} are not interlaced as 
in (3). 

We start by introducing some new notation that will be used throughout 
this section. If {xj} and {yj} is a solution of (5) and {p, q) is the corresponding 
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polynomial solution as defined in (13, 14), we can introduce the new polynomials 
Pr{z) = z^^'pil/z) and qr{z) = z"-«q{l/z) to describe the solution. Defining 
them by continuity at ^; = 0, we have 

Pr{z) = {z - Xi) ■ ■ -{z - XnJ, qr{z) = {z - yi) ■ ■ ■ {z - ynj ■ (39) 

Furthermore, we assume that the number of distinct roots of Pr (a;j-branch 
values) is n. We also order the roots such that we can write 

Pr{z) = {Z- Xi)l+''^ {Z - X2)l+''^ ■■■{Z- Xnf+'^- , 

where 1 + rij is the multiphcity of the root Xj , so that 

n 

nx = Deg(pr) = n + ^'qt. 

We start the analyis with a Lemma giving explicit expressions for the ak 
values. 

Lemma 4 For k >0, 

Proof: This result follows from an application of the residue theorem in complex 
analysis as follows. Let Cr be the circle in the complex plane with radius r. Since 
the roots oi p{z) are non-zero, the fmiction q/p is analytic within and on if 
£ is taken small enough, and the Cauchy integral formula gives 



ak 



l_ d'' qjz) 



, k>0, 1 r q{z) ^ 

z=0 =—4> , '/k+l '^Z. 



Setting 



0, fc < 0, 27ri Jc^ p{z)z 



f(.^ ._ ^r{z) _ z^y ^^q{l/z) 

^^'^■~Pr{z)- P{l/Z) ■ ^^^> 



and changing variable z ^ 1/z we get 

= ^ f^^ = ^ -^dz = — f^^^^ z f{z 

Hence, a„j^_„^+i+fe is given by the sum of the residues of f{z) (assuming we 
take small enough e). By (41) and the restriction /c > we see that its poles 
are located at the a;j-values and they have multiplicities 1 + rjj at Xj . Then 
(40) follows from the residue formula for a pole of a function g{z) at z* with 
multiplicity + 1 , 

1 d'^ 

Res(9,^*) = - hm^ —{z - z*f+^g{z). 

TJl z^z* dZ^ 

□ 
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When the branch values {xj } are distinct the expression for the Ofe elements 
simplifies. They can then be expressed as sums of the powers of {xj} in a way 
similar to the moments rrifc, but with weights different from one. We can also give 
a more concise description of the matrices Aq and Ai, which can be factorized 
into a product of Vandermonde and diagonal matrices. More precisely, we let 
V be the Vandermonde matrix 





( 1 


1 


1 \ 


v = 


Xi 

xl 


X2 













and introduce the diagonal matrices, 
W = 



V 



X = 



/Xi 



where Wj are the weights defined as 



w 



(42) 



(Note that Pr has only simple roots when {xj} are distinct, so Pr{xj) ^ 0.) 
Then we can show 

Proposition 4 If {xj} are distinct, then for k >0, 



(^ny—Ux+l+k — ^ ] 



and 



AiR = VWY"^, AoR = VWXV'^. 



(43) 
(44) 



Proof: When {xj} are distinct r]j = for all j and the expression (40) for the 
Sj-residue simplifies. 



lim 



{z - Xj)z^qr{z) ^ x'jqrjxj) 

Pr{z) p'^ixj) 

This shows (43). For (44) we set bk = a„ „^+i+fc. Then 



Ai—rR = 



( br 
br+1 



br+1 
br+2 



\br+nx W+Tix + l 

Prom (43) we then have, for fc > 0, 

f bk \ 
bk+i 



W+Tlx + l 



br+2nx I 



r = 0,l. 



V>k+nx/ 



Tlx 

^^E'^j'^ 


X] 


= V 


1 Wl4 \ 

W2X2 


= VW 


X2 








\WnxX^nJ 
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Consequently, 



Ai-rR = VW 



r+l 



r+l 



which concludes the proof. □ 

We now consider the implications of a positive definite ^i-R. It turns out 
that this is a necessary and sufficient condition to guarantee both distinct {xj} 
values and positive weights. We get 

Theorem 4 The matrix AiR is symmetric positive definite if and only if {xj} 
are distinct and the weights are strictly positive, Wj > for j = 1, . . . ,nx- 

Proof: We use the same notation as in Lemma 4 and set 



1 



Sjiz) = —{z-Xj) 
Vj- 



1+% 



Pr{z)' 



Wc note that Sj{z) is smooth and regular close to ^; = xj. Then by Lemma 4, 
for fc > 0, 



a. 



Next, we let v = {vi, . . . ,Vn^)'^ be an arbitrary vector in and recall that 
Pv{z) is the corresponding — 1 degree polynomial 



Py{z) =Vi+V2Z-\ h Vn^Z 



rix — l 



Then 



j=l k=l 



j=i fe=i e=i 



y lim ^S({z)yy z^+''-\vk = y hm ^Si{z)Py{zf 
'f^z^xidz'^i ^ '^f^ f^^^a;, d^w «w «v ; 

e=i j=i fe=i e=i 



If 



we can take 

P^z) = (z - {z - X2)i+^= • • • (z - 

Since 2(1 + ry^) = 2 + 2L%/2J > 2 + 2(%/2 - 1) > r]e and 



(45) 
(46) 



= 0, < ^ < fc. 
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for all smooth enough /(^), we get v^AiRv = 0, which contradicts the positivity 
of AiR. Hence, 



n 



n 



n + '^[r]j/2\ > Tlx - 1 = n + y^^T]i - 1. 
j=i e=i 



Since for any integer n > we have [n/2j < n — 1 it follows that all r]( = and 
h = rix- Hence, if AiR is positive definite, then {xj} arc distinct. 

To show the theorem it is now enough to show that, when {xj} are distinct, 
Alii is positive if and only if the weights are positive. Prom (45) we then have 



Clearly, when all we > 0, this expression is positive for v 0, and AiR is 
positive definite. To show the converse, we take P„ {z) to be the Lagrange basis 
polynomials Lj{z) of degree rix — 1 defined as 



This can be done for each j, which concludes the proof. □ 

We can now relate our conclusions with those in Markov's Theorem 3. We 
consider all solutions to (5), instead of those given by the integral relation (37) 
with a piecewise continuous function f{x). The extra condition (38) is then 
automatically satisfied, and we note that the positivity of AiR guarantees a 
unique solution also in our space of density functions (6). We view this as a 
corollary of Theorems 2 and 4. 

Corollary 1 If there exists a solution to (5), then the matrix in (38) is singular. 
When Tlx — 

there is a unique solution to (5) of the form (3) if and only if 
AiR is symmetric positive definite. 

Proof: We start by proving the singularity of (38). By (ii) in Proposition 2 

a coefficient solution c = (cq, . . . , c„^)'^ = {cq^c^ Y' satisfies Ac = 0. Since 
A = (ao Ai) it remains to prove that cook+i + OqC = 0. This was already 
proved in (36). 

Next, we prove the "if" part of the second statement. If AiR is symmetric 
positive definite it is non-singular and by (i), (iii) and (v) in Theorem 2 the 
minimal degree solution exists and is unique and Xj ^ yi for all (If Xj = 
for some j, then there is no zero r/i-valuc since T)cg{q*) = n by point (iii).) By 
Theorem 4 the corresponding branch values {xj} are distinct. It remains to 
show that, upon some reordering, the {xj} and {yj} are interlaced as in (3). 

Order the a;j-valucs in an increasing sequence and let m,k be the number 
of y^-values such that yj < Xk- Clearly, is increasing and < < Uy. 



v'^AiRv = ^5^(a;^)P„(x^)2 = Y,wePv{xtf . 





If Aii? is positive then 
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Moreover, sgn{qr{xk)) = (—1)"" ™'= and since lim2^ooPr(^) > 0, we also have 
sgn.{p'j.{xk)) = (—1)"^"'^. Hence, by also using the fact that Uy = Hx, 

Sgn(Wfe) = (^-l)ny-mk+n.-k ^ (_i)m,+fc_ 

We conclude that nik + k is even, which implies that m/j is in fact strictly 
increasing. Then, for A; = 1, . . . , — 1, we have mk+i > ruk + 1 and 

f^x > > ruk +nx — k => ruk < k. 

Similarly, ruk >mi + fc — l>fc — 1, sofc — 1< ruk < k, and therefore 

2k-l<mk + k<2k. 

Finally, since nik + A: is even we must have rUk = k, which implies that the 
values arc interlaced. 

Wo now consider the "only if" part. If there is a solution of the form (3), 
then the {.Tj}-valucs arc obviously distinct and rrik = k. By Proposition 4 the 
weights are then given by (42) and they are positive since, as above, sgn(z/;fc) = 
(— 1)'"'=+'^ = 1. It follows from Theorem 4 that A\R is positive definite. □ 



7 Outlook 

Several interesting issues may be worth mentioning: 

1. Computational complexity in a finite difference implementation: one can 
consult the article [14] where practical implementation issues and several 
examples of increasing complexity have been addressed in the context of 
geometric optics problems. In particular, comparisons with Lagrangian 
(ray-tracing) solutions are shown. 

2. Extension to higher dimensions: nothing seems to be done in this direction 
at the time being; see however the last sections of [20] and the routines 
based on complex variables in [11, 9] for "shape from moments". 

3. A very special case of the trigonometric moment problem can be solved 
by means of a slight variation of the algorithms presented here, in [14] and 
in Section IV. A of [9]. That is to say, one tries to invert the following set 
of equations: 

n 

jjLj exp{ikXj) = nik, k = 0,...,n. (47) 

Let us state that in case the n + 1 real frequencies Xj are known, the set 
of complex amplitudes fxj are found by solving a Vandermonde system: 



/ 1 


1 


\ 


/ Mo 


\ 


f 


mo 


\ 


exp(iAo) 


exp(iA„) 








mi 


\ exp(mAo) 


exp(mA„) 


/ 




J 




m„ 
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The frequencies can be found through a byproduct of [9, 14] as we state 
now: let us suppose n is odd {i.e. the number of equations is even), we 
form the two matrices, 

/ mi ••• rrin+i \ 

2 

\ nirz+i ■ ■ ■ rUn I 

and then the frequencies can be obtained through a generahzed eigenvalue 
problem, AiVj = XjA2Vj, j = 0, ...,n. This kind of algorithm can be 
used to check the accuracy of the classical FFT and will be studied in a 
forthcoming article. 
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